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Using a variational numerical method we compute the time-evolution of the Hubbard-Holstein 
bipolaron from its ground state when aX t = Q the constant electric field is switched on. The system 
is evolved taking into account full quantum effects until it reaches a quasistationary state. In the 
zero-field limit the current shows Bloch oscillations characteristic for the adiabatic regime where 
the electric field causes the bipolaron to evolve along the quasiparticle band. Bipolaron remains 
bound and the net current remains zero in this regime. At larger electric fields the system enters the 
dissipative regime with a finite quasistationary current. Concomitantly, the bipolaron dissociates 
into two separate polarons. By examining different parameter regimes we show that the appearance 
of a finite quasistationary current is inevitably followed by the dissociation of the bipolaron. 
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I. INTRODUCTION 



Nonequilibrium phenomena in interacting many-body 
systems have recently drawn significant attention. One of 
the most intriguing challenges is to understand how a sys- 
tem responds to an external electric field and to unravel 
the microscopic mechanism that gives rise to a constant 
current and a steady increase of the total energy of the 
system. From a theoretical point of view, a non-linear 
current-voltage characteristics with a threshold behav- 
ior is well-establishedfi"— and recently the appearance of 
a regime of negative differential resistivity has been ob- 
served in several interacting systems Moreover, the 
onset of Bloch oscillations (well-known for the nonin- 
teracting systems) has been observed in different corre- 
lated systems at very large electric ficlds^i^i^ and even 
in the integrable modelsi^ Oscillations in the current re- 
sponse can also exhibit a more complex pattern leading 
to beats in a certain regime of model parameters i^iii Nev- 
ertheless, the heating arising from the energy flow to the 
system represents a serious obstacle in calculation of a 
constant nonzero current for systems at half-filling ^liS 
therefore most of the calculations investigating current- 
voltage characteristics have been limited to one dimen- 
sionj^ii^ One of the widely used concepts is to couple the 
system to a thermal bath, which at a loss of quantum co- 
herence enables calculations of a quasistationary current. 
Recently, the influence of coupling to a thermal bath on 
the current-voltage characteristics of the Hubbard model 
has been invcstigatedii The reported influence is similar 
to the role of quantum phonons on dissipation of the 
excess energy and consequently on the current-voltage 
characteristics of the Holstein modelJ^ 

A complementary approach to investigation of quan- 
tum systems out of equilibrium is to study dynamics of a 
single carrier propagating in dissipative medium^i^ii^ in 
particular a charge carrier driven by a constant electric- 
fieldji^ Such problems have been recently solved for a 
driven carrier doped into Mott insulator— and a driven 



carrier coupled to Holstein phonons where the nonzero 
current arises due to constant emission of magnons and 
phonons, respectively. Recently, using a state-of-the-art 
truncated Lanczos methodi^""— the quasistationary cur- 
rent in a doped 2D Mott insulator coupled to phonons 
was calculatedfi^ The advantage of such approaches is 
refiected in calculation of a stable well-defined current 
while preserving the full quantum nature of the prob- 
lem. In addition, the total energy gained by hopping of a 
charge carrier along the direction of the electric field is en- 
tirely absorbed within the microscopic model. This idea 
enables the formulation of the time-evolution of different 
subsystems without coupling to an external reservoir. 

The aim of our study is to extend a recent study of a 
driven charge carrier coupled to Holstein phonons^ to a 
two-particle problem as described within the Hubbard- 
Holstein Hamiltonian. There are two sources of the 
particle-particle interaction in this problem: a) an indi- 
rect interaction, mediated by the electron-phonon cou- 
pling and b) the direct on-site Coulomb interaction. This 
motivates us to address the following questions: (i) Is it 
possible to drive a bound pair of electrons out of equilib- 
rium in such way that a finite current is obtained without 
breaking the pair? We show that a Hubbard-Holstein 
bipolaron at initial time t = 0, driven by a constant 
electric field, begins to dissociate as soon as a nonzero 
quasistationary current is reached, (ii) How is the qua- 
sistationary current calculated per particle changed in a 
two-electron system with respect to the single-electron 
case? (in) We analyze beats in the transient current re- 
sponse. We show that the beats in our model emerge 
if the energy scale associated with the Bloch oscillations 
within the quasiparticle (QP) band coincides with the 
energy gap between the QP band and the low-energy 
continuum of excited states. 

The paper is organized as follows. We introduce the 
model and numerical method in Sec. II. In Sec. HI we 
show numerical results for the short-time behavior of a 
driven Hubbard-Holstein bipolaron with a focus on the 
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dissociation of a bound state due to the constant external 
electric field. In Sec. IV we address properties of the 
current-field characteristics in the quasistationary state 
while in Sec. V we analyze the Fourier spectrum of the 
real-time current. We give conclusions in Sec. VI. 



II. MODEL AND NUMERICAL METHOD 

We define the one-dimensional time-dependent 
Hubbard-Holstein Hamiltonian, threaded by an external 
flux 

La j 

j j 

where and are local electron and phonon creation 

operators on site j, respectively, and ??j = J2cr a'^j<^ 
is electron density, ujq denotes the dispersionless opti- 
cal phonon frequency, to is the nearest-neighbor hopping 
amplitude and U represent Hubbard repulsion between 
on-site electrons. We set wq/^o = 1-0 throughout the pa- 
per. The dimensionless electron-phonon (EP) coupling 
strength is defined as A = g'^/2tQUJo. The electric field is 
applied through a time-dependent Aharonov-Bohm flux 
<j>{t) = —Ft, which generates force on the charged parti- 
cle due to the Faraday's law. The constant electric held 
F is switched on at time t = and is measured in units 
of [to/eoa], where cq is the unit charge and a is lattice 
constant. The time is furthermore measured in units of 
[to/^], and we set to = h = cq = a ~ 1. 

In our calculations we obtain reliable numerical results 
for electric fields in the range between 10^ — IQ^V/cm, 
if the model parameters are set to to — O.leF and 
a = 10^^*^771. These fields become relevant for compar- 
ison with experiments as presented in Ref.— if the ex- 
trapolation to zero temperature T ^> is taken. 

We solve the time-dependent Schrodinger equation for 
two electrons coupled to phonon degrees of freedom, 
with the use of the numerical method based on the 
exact diagonalization of the variational Hilbert space 
(VHS) that led to numerically exact solutions of the po- 
laron and bipolaron ground and low-lying excited-state 
properties<^2i2ii2i In the case of the Hubbard-Holstein 
model the method generates the VHS starting with the 
translationally invariant initial state \ipo) where both 
electrons are located on the same site with no phonon 
degrees of freedom. The VHS is then generated by re- 
peatedly applying the off-diagonal terms of Hamiltonian 
in Eq. ©, 

{i^f-^^)} = [ij,„ + E(^.n^'i¥'o), (2) 

m— 1 

where Hkin and Hg are the first and the second term of 
the Hamiltonian in Eq. ([T]). Parameters Nh and M deter- 
mine the size of the VHS. The parameter Af > 1 ensures 



good convergence in the strong EP coupling regime. To 
reach the intermediate or weak coupling regime A < 1 , we 
introduced an additional parameter Nphmax which limits 
the maximal number of phonon quanta, enabling larger 
values of Nh ■ 

We first calculate the ground state at F = of the 
Hubbard-Holstein Hamiltonian in Eq. ([T|). Equilibrium 
properties of bipolarons have been extensively studied 
in the literature^"— with particular emphasis on the 
conditions for formation of a bound bipolaron state. We 
switch on the uniform electric field F at t = and start 
the time propagation using the time-dependent Lanczos 
technique^i The time step is taken to be small enough 
{dt/tB = 0.5 X 10"^ where ts = 2tt/F denotes the Bloch 
oscillation period) to ensure good numerical accuracy. 
We investigate only the Sf^j = subspace. An important 
observable is the time-dependent average of the current 
operator per particle, j{t) = {I{t))/2, where 



lit) 



-iFt J 
^1 , 



H.c. 



(3) 



Since we are dealing with the time-independent field F, 
the time integral of the current is directly related to a 
change of the total energy 



j{t')dt' = Ahit)/2F, 



(4) 



where Ah{t) = {H{t)) - {H{t = 0)),^ We also calcu- 
late the average number of phonon quanta at time t, 

{nph){t) = (Ej 

Due to the finite size of the VHS we are able to 
run the propagation only until the system reaches 
the boundary of the variational space. There are two 
typical boundaries the system can reach: (i) the size 
of the bipolaron becomes larger than iV^, or (ii) the 
total number of phonons in the system approaches 
the maximal allowed value Nphmax- We checked the 
finite size effects by comparing the expectation values 
of different observables for different system's sizes. In 
Fig. [1] we represent the time dependence of the relative 
difference between energies for different systems sizes, 
namely \Ei — -EjVh=i5|/-E'Wh=i5' where i represents the 
parameter Nh of the functional generator in Eq. ^ 
and the maximum number of phonons were limited to 
Nphmax = 15. Our simulations were stopped when the 
relative energy difference between largest and second 
largest system was more that 5%. Using the sum rule, 
Eq. we were able to independently control the 

precision of the time-evolution, where errors occur due 
to time discretization. 



III. 



REAL-TIME BEHAVIOR 



Previous studiesii^i^ have shown the importance of 
the low-energy spectrum for the short-time behavior of 
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FIG. 1: (Color online) Time dependence of the relative dif- 
ference between energies for different systems sizes |i?jv^ — 
-EislZ-Bis, where Bjv^ is energy of the system generated by the 
functional generator, see Eq. ((2]), using Nh as parameter and 
maximal number of phonons were limited to Nphmax = 15. 
We set A = 0.9 and F = 1.0 in this figure. 

the interacting systems after switching on the static elec- 
tric field. In terms of short-time dynamics there exist 
at least three distinct energy scales. First is the bind- 
ing energy of the bipolaron i5 representing the energy re- 
quired to dissociate a bound bipolaron into two separate 
polarons, S = Eu — "^Epoi, where E^i and Epoi are the 
energies of a bipolaron and a polaron, respectively. The 
second is the phonon frequency wq, the third is the gap 
in the excitation spectrum A separating the top of the 
bipolaron band from either the continuum of states or the 
excited bipolaron band. Due to the existence of a finite 
gap, there exists a threshold field Fth separating the adi- 
abatic regime for F < Fth and a dissipativc regime with 
a finite current at F > Fth, see also Refsjii^. Character- 
istic features of the adiabatic regime are periodic Bloch 
oscillations as the bipolaron moves along the bipolaron 
band. As a consequence, a time-averaged current in the 
adiabatic regime is zero. On the other hand, in the dis- 
sipativc regime a quasistationary state is reached after 
a short transient time t > ttr- A quasistationary state 
is characterized by a constant time-independent current 
= J and a linear growth of the total energy of the 
system, 

Ah{t)^j2F. (5) 

In the case of a driven polaron it has been shown that the 
steady increase of Ahit) is entirely due to the increase 
of phonon excitations left in the wake of the traveling 
polaroni^ 

A reasonable conjecture in the case of the bipolaron 
is that for a large binding energy, i.e., for S > ujq, the 
system on short time scale propagates as a bound bipo- 
laron, depositing the excess energy in the form of excited 
phonons left in its wake. Our calculations show that in 
the weak coupling regime (A <C 1) the bipolaron dissoci- 
ates on a time scale much shorter than is . We therefore 



focus in the following on the intermediate and strong 
coupling regimes (A = 0.5 and 0.9, respectively), where 
the propagation of a bound pair under the electric field is 
more likely. However, as we shall show, such propagation 
is not observed in our calculations, instead, the bipolaron 
begins to dissociate as soon as the system enters the dis- 
sipativc regime with a finite electric current. 

In this section we focus on U = 0. The motivation for 
choosing A = 0.5 and A = 0.9 is that the binding energy 
is smaller than ujq {5 sa 0.58 < ujq) for A = 0.5 while 
for A = 0.9 it is much larger than ujq {S « 2.04 > ljq). 
We also note that in the case of A = 0.9 there exists 
an excited bound state at A ^ 0.85 above the ground 
state at fc = 0. This state represents a bound state of 
the bipolaron with excited phonon cloud. Similar bound 
excited states have been well-studied for the case of the 
Holstein polaron.i^^^ They form dispersive bands just 
below the one-phonon continuum that starts at luq above 
the ground state energy. 

For A = 0.5 and electric field F ^ 1/7 the system dis- 
plays nearly adiabatic Bloch oscillations in the real-time 
current j{t) as shown in Fig. [^IJa), which are consistent 
with the bipolaron ground state dispersion. Weak damp- 
ing due to inelastic scattering on phonons is observed, 
that is in turn reflected in slow increase in the total sys- 
tem energy Ah{t) in Fig. [DJc). For larger electric field 

= 1/3 the reminiscence of the Bloch oscillations are 
still seen, however the current remains positive at all 
t > 0. From the available data we estimate the thresh- 
old field Fth ~ 0.2. Although for A = 0.5 we have not 
strictly reached the quasistationary state since oscilla- 
tions in current remain well pronounced, the total energy 
Ah{t) shows increase in time which is approximately lin- 
ear. The increase of Ah{t) allows us to use a linear fit 
according to Eq. (O and to calculate the value of qua- 
sistationary current j. The corresponding current-field 
characteristics and the influence of model parameters on 
J will be discussed in Sec. IIVI 

We also observe a quasi-linear increase of the total 
phonon quanta in the system {nph){t), see Fig.[2jb), and 
we find that Ah{t) > Ldod{nph) /dt. While the equality 
would mean that the entire energy absorbed from the 
electric field is absorbed by the lattice, the inequality 
signifies that there exists excess energy that is responsible 
for the dissociation of the bipolaron. Further insight into 
this process can be obtained by calculating the average 
particle distance d{t) defined as 



d{t) = (6) 

y ij 

shown in Fig. El^d). We detect an overall increase of d{t) 
for both values of F, signaling the dissociation of the 
bipolaron. We discuss the process of dissociation in more 
detafl in Sec. ImXland Sec. HITbI 

In the strong-coupling regime at A = 0.9 we calculate 
the real-time response of the system at F = 0.5 « Fth, 
shown in Figs. [2Ia), EJc), EJe) and [21 g)- Remarkably, 



FIG. 2: (Color online) j{t), {nph){t), Ah{t) and d{t) vs t/fs 
for A = 0.5. We represent characteristic fields in nearly adia- 
batic regime and when F > Fth, i.e., F = 1/7 and F = 1/3. 
The accuracy of the propagation was checked by the compari- 
son to the energy-gain sum rule, Eq. (U). Parameters defining 
the functional generator of Eq. ((2| for A = 0.5 used through- 
out the work were A^^ = 17, M = 2, Nphmax ~ 15. 



beside the previously mentioned effects, we observe beats 
on a longer time-scale. We discuss their origin in Sec.lVl 
If the values of electric field arc increased above F « 
Fth, the beats become less pronounced. For F = 1 and 
F ~ 2 wc obtain, after initial oscillations, a well defined 
quasistationary current j(t) and a nearly linear time- 
dependence of {nph){t) and Ah{t) characteristic for the 
dissipative regime, see Figs. [31Jb), [31[d) and[31Jf). 



FIG. 3: (Color online) j{t), {nph){t), Ah{t) and d{t) vs t/ts 
for A = 0.9. Left column: F = 0.5 and a maximal t/tB ~ 15. 
Right column: F — 0.5, 1 and 2 with a maximal t/ts ~ 8. 
Thin horizontal line in (b) indicate the quasistationary cur- 
rent J. Thin lines in (e) and (f) represent a linear depen- 
dence of Ah{t), which indicates the quasistationary state. 
Parameters defining the functional generator of Eq. for 
A = 0.9 used throughout the work were Nh = 13, M = 3 and 

Nphmax 15. 



ited with the maximal allowed size of the Hilbert space, 
our results strongly suggest that in the limit when t ^ oo 
the bipolaron dissociates into two separate polarons. 

Additional information about the bipolaron dissocia- 
tion can be obtained if d{t) is compared to the average 
distance traveled by the center of mass of two particles 
Ax{t), where 



Ax{t) = 



Ah{t) 
2F 



(7) 



A. Bipolaron dissociation 



Using Ax{t) it is convenient to define a ratio 



The central goal of our work is to understand whether 
in the long-time limit the bipolaron remains bound (in 
some parameters regimes) as it travels under the influ- 
ence of the constant electric field, or it dissociates into 
two separate propagating polarons. The most interesting 
is the regime of strong EP interaction where the bipo- 
laron binding energy 6 is larger than the phonon energy 
wq. Such case is shown in Figs.[3lg) and[3ljh) for A = 0.9. 
The average distance d(t) is steadily increasing even in 
the near-adiabatic regime, and the slope of the linear in- 
crease of d(t) is enhanced in the dissipative regime, as 
seen for F = 1 and 2 in Fig. [3jh). 

We should stress that the overall increase of d{t = 
^max) after the maximal propagation time relative to the 
initial ground state distance d{t = 0) is more than 6-fold 
in both cases, i.e., at F = 1 and F = 2. Even though the 
maximal propagation time is in our time-evolution lim- 



Ax{t) 
Ad{t) ' 



(8) 



which represents the ratio between the average distance 
Ax{t) traveled by two particles along the field direction, 
and the relative increase of the distance d{t) between the 
particles defined as Ad{t) = d{t) — d{t = 0). 

If the bipolaron propagated without dissociating into 
separate polarons, we would expect a linear increase of 
ri{t) in the quasistationary state. Instead, results for 
both A = 0.5 and 0.9 displayed in Fig. |4] show a clear 
tendency toward a constant value of the order of one, in- 
dependently on the strength of F. For small F ~ Fth, 
T]{t) shows damped Bloch oscillations while for larger F 
the approach toward a constant is rather monotonous. 
This results suggest that the dissociation of the bipo- 
laron and the appearance of the quasistationary current 
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FIG. 4: The ratio ri{t) as defined in Eq. ©, vs t/tij. Up- 
per panel and lower panel correspond to A = 0.5 and 0.9, 
respectively. 



emerge simultaneously as the system evolves from a tran- 
sient regime to the quasistationary state. 



B. Phonon correlation function 

We also compute the average number of phonon quanta 
located at a given distance r from (both) electrons. To 
study this effect we define 



7(r) = 



1 



(9) 



fulfilling the sum rule ^^^jir) = 1. Correlations func- 
tion 7(r) for A = 0.5 (0.9) and F = 1/3 (1) arc displayed 
in Fig. [5l giving additional insight into the dynamics of 
the moving bipolaron. The most characteristic feature 
is a pronounced asymmetry of 7(r) with respect to the 
electron positions at r = that grows with time. The 
asymmetry emerges due to phonon excitations extending 
behind the moving particles, which contain the excess 
energy absorbed from the external electric field. Here we 



FIG. 5: (Color online) Zo.gio(7(r)) for A = 0.5 (0.9) and 
F = 1/3 (1) computed at different times. Note that the log- 
arithm of the correlation functions is presented. The vertical 
dashed lines represent the average distance between the elec- 
trons d at a certain time. The arrow represents the direction 
of the moving bipolaron. Green dashed-dotted lines represent 
the corresponding logio{'y{r)) of the Holstein polaron model.— 
We show logio{'y{r)) for \r\ ^ A^^ where results for different 
system sizes show good convergency. 



note that the electron is moving in the direction of r > 0. 
This effect is similar to the polaron case.— The second fea- 
ture is the increased amount of phonon excitations in the 
forward direction, which arises due to two distinct con- 
tributions. While the first contribution is due to damped 
Bloch oscillations analogous to the polaron propagation^^ 
the second contribution is entirely due to a two-particle 
effect, in which case there is always one particle travel- 
ing ahead in the electric field. It is the second particle 
which follows the first one that in turn detects phonon 
excitations left by the first particle, thus generating extra 
weight of 7(r) in the region r > 0. Therefore, to obtain 
the complete explanation of the increasing 7(r) in the 
forward direction one has to take also into account the 
increasing average distance between two particles. 

For comparison we added results for the polaron case, 
presented by dashed-dotted lines in Fig. [SJ which show 
a farther extend of 7(r) behind the moving polaron and 
a smaller reach of 7(r) in the forward direction, the lat- 
ter being consistent with the argument given above. A 
shorter phonon disturbance behind the traveling bipo- 
laron indicates a smaller center of mass velocity in com- 
parison to the polaron. We further elaborate on this issue 
in the following section. 



IV. QUASISTATIONARY CURRENT 

We next present results for the quasistationary current 
J, its dependence on F and U, and comparison with the 
polaron case. In Fig. HJa) we display the current-field 
characteristics for A = 0.5. We have limited our calcula- 
tions to commensurate values of = uio/n for F < ujq 
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and F = uluq for F > ujq and n is integer. The cur- 
rent in the bipolaron case is in this case smaller than 
the current in the polaron case. This is a consequence of 
the fact that in the former case one of the two electrons 
is traveling in the wake of the other electron where the 
wake consists of the finite number of phonon excitations. 
The intuitive explanation is that in this region the ef- 
fective temperature is elevated that in turn leads to the 
lowering of the current. To check this assumption we per- 
formed a numerical test, where we have propagated the 
polaron from the ground state until it has reached the 
quasistationary state. Then we changed the direction of 
the electric field, which prompted the polaron to propa- 
gate backwards into the phonon rich area. After initial 
oscillations the polaron has reached the quasistationary 
state with a lower current, i.e., for A = 0.5 the current 
in the phonon rich area was roughly half of the value, 
obtained when it traveled into the undisturbed region. 

Assuming that j is reduced in a two-particle case due 
to the arguments given above, and setting r]{t) w 1 in 
the quasistationary state for A = 0.5 (see Fig. UJ^a)), it is 
possible to estimate the j-F characteristics of a driven 
bipolaron from a polaron j-F characteristics. If Ax{t) = 
-|- a)/2, where vi is the polaron velocity in the qua- 
sistationary state and a < 1 is the renormalization of the 
second particle's velocity, then Ad{t) = vi{l — a). Using 
a rough estimate 77 w 1 yields the renormalization factor 
a w 1/3. This enables us to estimate the two-particle 
current as ]2{F) ~ Indeed, these results shown 

in Fig.in^a) are in reasonable agreement for F ^ 0.5 with 
the quasistationary current of a driven bipolaron calcu- 
lated from the real-time dynamics. Using these consid- 
erations we can understand almost linear dependence of 
d{t) in Fig. [21[d) and[31^h), since the particle the follows 
the first one moves with a renormalized but a constant 
velocity. 

At larger lambda A = 0.9 as shown in Fig. H^b), the 
bipolaron remains in the near adiabatic regime up to the 
threshold field Fth ~ 0.5 that is much larger than Fth ^ 
0.1 for the polaron case. The dependence of Fth on model 
parameters is better understood when loosely applying a 
simple Landau-Zener (LZ) formalis m'^^i'^^ to a multi-band 
energy spectrum of the bipolaron, which givesii^ 



pLZ 



(A/2)2 



(10) 



where v is the group velocity of the quasiparticle. Much 
larger Fth in the case of bipolaron is attributed to 
the much narrower bipolaron bandwidth Wbi in com- 
parison to the polaron Wpoi, which leads to a larger 
gap A and a lower group velocity v of the quasipar- 
ticle. We can estimate the scaling of F^^ for bipo- 
laron in the strong coupling limit. The gap is A « 
a;o(l — 4<:oexp(— 8Ato/wo)) and v is approximated as max- 
imal value of dEfii/dk, which in the asymptotic expansion 
scales as V K, 2^0 exp(— 8Ato/wo). Therefore to the lead- 
ing order in A, Fth scales as Fth ~ (i^o/^o) cxp(8Ato/i^o)- 
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FIG. 6: (Color online) (a) Quasistationary current j vs elec- 
tric field F at A = 0.5 for a polaron (Ih) and a bipolaron 
(2h). The bipolaron quasistationary current estimated from 
the polaron quasistationary current is denoted by j2 (see text 
for details). The uncertainty of J is represented with error 
bars and is a consequence of transient region; see discussion 
in text, (b) Quasistationary current J vs electric field F for 
A = 0.9 and different values of the Hubbard repulsion [/ = 0, 
(7=1, and (7 = 2. The uncertainty for J is represented 
with error bars and is again a consequence of transient re- 
gion, (c) The threshold electric field Fth as a function of the 
electron-phonon coupling A for (7 = 0. The inset represents 
dependence of the threshold electric field Fth as a function of 
electron-electron repulsion U at fixed A = 0.9. 



To extract Fth from ] we used 



(11) 



where cto and Fth are fitting parameters. 

Nearly linear scaling of log(Ft;i) with A is represented 
in Fig.lHKc). Here we should emphasize that in the leading 
order of A, Fth is increased predominantly due to the 
narrowing of the bandwidth and consequently reducing 
of the group velocity v in Eq. ([TU)) . 

In Fig. mb) we also present the influence of the Hub- 
bard repulsion \] on the quasistationary current. The 
most important effect is the increase of j with increas- 
ing V for F < 2. In several recent studies of interacting 
systems driven out of equilibrium, an increase of a qua- 
sistationary current due to stronger correlations has been 
observed i^ii^ In the case of the Hubbard-Holstein bipo- 
laron the effect of V is to reduce the double occupancy 
and consequently the binding energy b. Moreover, in the 
framework of the strong EP coupling limit, V reduces 
the effective mass and thus increases the bandwidth)^ 
which in turn leads to a decrease of the gap in the bipo- 
laron excitation spectrum. These combined effects lead 
to a predominantly linear decrease of Fth with increas- 
ing IJ in the regime U ^ 1, as represented in the inset 
of Fig. [HKc). Nevertheless, Fth tends to a constant value 
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Fth ~ 0.2 already around U ^ 1.5 while the bipolaron in 
the ground state remains bound up to Uc ^ 3, see Ref»^. 
We also note that a similar threshold value Fth ^ 0.1 — 0.2 
is found in the case of a single polaroni^ The bipolaron 
thus reaches threshold field values of a single polaron even 
in the regime when U < Uc- 

Here we should note that in case of A = 0.9 we were 
unable to reach (for all electric fields) a fully dissociated 
state of two separated polarons due to a long dissocia- 
tion time and a limited variational Hilbert space. Even 
in this transient regime system shows clear characteristic 
features of quasistationary dependence, i.e., linear de- 
pendence of energy with time, see Fig. [3jf). j-F depen- 
dence for A = 0.9 represented in Fig. [HJb) is calculated 
in this regime, where the effects of the correlations be- 
tween particles are still strong and we would expect that 
after the transient regime the system would enter a new 
regime, when bipolaron would be well dissociated and 
dynamics would be more similar to the A = 0.5 case. 



V. FOURIER ANALYSIS 

Following the idea of Refsj^*^ we perform the Fourier 
analysis of j(t) in order to extract additional information 
about the dynamics of the driven system. In doing so we 
first subtract the linear and quadratic terms from jit) 
data in all cases where applicable. In Fig. [7] we repre- 
sent the normalized Fourier transform (the integral over 
whole spectrum is 1) of the real-time current for A = 0.9 
and different electric fields. We are measuring time in 
units of Bloch time 2-k /F and the angular frequency is 
accordingly given by w = a; * is . 

In the (nearly) adiabatic regime, for F = 1/3, the bipo- 
laron exhibits Bloch oscillations, with the doubled Bloch 
frequency W2b = 2Cjb — 47r. This is refiected in a well 
pronounced peak in the Fourier spectrum of i{t) at Cj2b- 
We as well observe well separated peaks at higher har- 
monics of a)2B, however, with decreasing amplitude. 

In the strong-coupling picture processes responsible 
for such peaks are diagonal transitions,— i.e., processes 
where bipolaron exhibits coherent propagation through 
the lattice. For example, if the bipolaron jumps as a 
bound pair for one site along the field direction, the en- 
ergy difference between the initial and the final state is 
2F. The related frequency for such a process is w = 2F 
or = 47r. Jumps for more sites represent higher order 
processes, which are responsible for generation of higher 
harmonics. Other peaks are related to different energy 
scales in the system, such as wq and A as indicated in 
Fig. [71 

For F = 1/2 we noticed pronounced beats in different 
observables, see Figs. I^a), (c), (e) and (g). In Fig. [7] at 
F = 1/2 the most pronounced peak corresponds again 
to w = Att, however, it this case the bipolaron Bloch fre- 
quency (jJ2B matches the phonon frequency ujq. This is in 
contrast to the case at F = 1/3 where we notice no such 
accidental degeneracies, and no beats appear in j{t). An- 
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FIG. 7: (Color online) Normalized Fourier transform of the 
real-time current j{t), where linear and quadratic terms have 
been extracted. We set electron-phonon coupling A = 0.9 
and electric fields F = 1/3, 1/2, 1. Arrows mark the angular 
frequencies corresponding to the distinct energy scales of the 
system. The horizontal dashed line represent the frequency 
corresponding to the bipolaron binding energy. 



other condition for the appearance of the beats is that the 
electric field is in the vicinity of the threshold field, i.e., 
F Fth, since at larger F increased energy dissipation 
to phonons overdamps Bloch oscillations. Pronounced 
beats in j{t) were observed as well by other authors ex- 
ploring different driven models.— 

For F = 1 the system reaches the quasistationary state 
in a time t ^ ts with a nearly linear total energy in- 
crease A/i(t), see Fig. [211), and a rapid increase of d{t), 
see Fig. [3Kh). For this reason the Fourier spectrum in 
Fig. [7] of j{t) shows less pronounced, broader peaks. We 
note that the frequency corresponding to the bipolaron 
binding energy 6 as marked with the dashed line is com- 
parable to U!2B- 

The influence of the Hubbard repulsion on the dynam- 
ics of j(t) is presented via the Fourier transform in Fig. [5] 
for the case of A = 0.9 and F — 1/2. The spectrum at 
finite U = 1.0 shows broader peaks near the character- 
istic frequencies in comparison to C/ = results. The 
effect of increasing U is in agreement with the consider- 
ations given is Sec. IIVI We again note that the spectrum 
is broadened when the frequency corresponding to S be- 
comes comparable to ll)2b- Despite much broader spec- 
trum a,t U = 1, J remains small, comparable to the U = 
case. At larger U = 2, the peaks in the spectrum become 
indistinguishable, while simultaneously a significant qua- 
sistationary current appears, as seen in Fig. IHI^b). 



VI. CONCLUSIONS 

We have studied the time-evolution of the Holstein- 
Hubbard bipolaron from its ground state at zero temper- 
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FIG. 8: (Color online) Normalized Fourier transform of the 
real-time current j{t), where linear and quadratic terms have 
been extracted, for different values of Hubbard repulsion U. 
We set A = 0.9 and F = 0.5. Arrows mark the angular 
frequencies corresponding to the distinct energy scales of the 
system. The horizontal dashed line represent the frequency 
corresponding to the bipolaron binding energy. 



ature after a constant electric field has been switched on 
at i = 0. Using an efficient variational method, defined 
on an infinite one-dimensional chain, we time-evolved the 
many-body Schrodinger equation while preserving the 
full quantum nature of the problem, until the system 
has reached a quasistationary state. In the limit of small 
electric field the bipolaron evolves along the quasiparticle 
band while the current shows characteristic Bloch oscil- 
lations. In this regime the propagation is adiabatic, there 
is no increase of the total energy, the net current remains 
zero and the bipolaron remains bound. At larger electric 
fields we detect an overall increase of the total energy, 
a finite net current appears and the system enters the 
dissipative regime. Due to dissipation of the gained po- 
tential energy into lattice vibrations, Bloch oscillations in 
this regime become damped. After a transient time the 
system enters a quasistationary state with a constant cur- 
rent. We note that there is no clear distinction between 
the adiabatic and dissipative regime. While a true adi- 
abatic regime exists only in the limit when F — > 0, a 
sizable net current signaling the dissipative regime ap- 
pears for F > Fth- 

The focal result of our work is that in the dissipa- 
tive regime the bipolaron dissociates into two separate 
polarons. By examining different parameter regimes we 
realized that the appearance of a finite quasistationary 
current is closely connected with the dissociation of the 
bipolaron. Even though our calculations were limited to 
large F and short propagation times, our results strongly 
suggest that bipolaron has a finite lifetime as long as 
the system displays a non-zero j. This hypothesis is 
supported by the calculation of rj^t) that in the large- 
time limit tends to a constant, largely independent of F. 



This result demonstrates that the dissociation and the 
appearance of the quasistationary current emerge simul- 
taneously as the system evolves from a transient regime 
to the quasistationary state. 

Here we should comment that due to dispersionless 
phonons the bipolaron can gain energy when two par- 
ticles occupy the same site, but to achieve a non-zero 
quasistationary current charged particles have to start 
hopping along the direction of electric field, and in doing 
so breaking the bound pair. In Ref4i the Green function 
study at finite temperature (without electric field) calcu- 
lated the bipolaron dissociation time, which goes to zero 
for dispersionless phonons, while this is not the case for 
phonons with dispersion. For the problem studied here 
this may indicate that phonons with dispersion would not 
always lead to dissociation of bipolaron as soon as a finite 
quasistationary current is obtained. Further research in 
this field is necessary to resolve the issue. 

A linear time-dependence of a real-space expansion 
of particle densities or related quantities has been ob- 
served in several recent nonequilibrium problems 
Even though we study a two-particle problem on an infi- 
nite lattice and we do not deal with a well-defined ther- 
malization, we can consider a cloud of emitted phonon 
excitations as a subsystem with an elevated effective tem- 
perature. The quantum Monte Carlo study^ of the ther- 
mal dissociation of the Hubbard-Holstein bipolaron pro- 
vides evidence that a bipolaron dissociates as the tem- 
perature becomes comparable to the binding energy S. 
In comparison, our results indicate that the bipolaron 
dissociates as soon as a finite net current appears around 
F Ffh- This occurs even in the case of 5 > wq, which 
is realized for A = 0.9 in our calculations, see Fig. [3t^h). 
We have also checked the time dependence of systems 
at even larger electron-phonon coupling, up to A < 1.5. 
With increasing A the qualitative behavior remains the 
same, the bipolaron as well dissociates with the disso- 
ciaton time that is increasing with A and the threshold 
field is also roughly exponentially increased as shown in 

Fig. lire). 

Yet another noteworthy finding reveals that the qua- 
sistationary current per particle is in case of the disso- 
ciated bipolaron smaller than in the single polaron case 
since the velocity of one particle is diminished due to the 
scattering on phonons generated by the motion of the 
other particle. Comparison of the bipolaron j-F charac- 
teristics at A = 0.9 with the polaron case (not presented) 
shows that a quasistationary current of bipolaron is not 
always lower than that of the polaron, which is possibly 
a consequence of remnant correlations between particles 
due to the incomplete dissociation of the bipolaron in the 
transient regime. 

To summarize, we have demonstrated the dissociation 
of the bipolaron under the influence of the electric field. 
This result predicts that as long as a given dilute sys- 
tem of bipolarons can be described by the short range 
Hubbard model in ID and the on-site coupling term be- 
tween the charge density and the dispersionless phonons. 
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bipolarons do not carry electric current, instead they dis- 
sociate into separate polarons as soon as the electric field 
is strong enough to yield a finite current. 
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